A Growth-based Optimization Algorithm for Lattice Heteropolymers 
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An improved version of the pruned-enriched-Rosenbluth method (PERM) is proposed and tested 
on finding lowest energy states in simple models of lattice heteropolymers. It is found to outperform 
not only the previous version of PERM, but also all other fully blind general purpose stochastic 
algorithms which have been employed on this problem. In many cases it found new lowest energy 
states missed in previous papers. Limitations are discussed. 



Lattice polymers have been studied intensively to un- 
derstand protein folding, one of the central problems of 
computational biology. A popular model used in these 
studies is the so-called HP model |l|, |2J where only two 
types of monomers, H (hydrophobic) and P (polar) ones, 
are considered. Hydrophobic monomers tend to avoid 
water which they can only by mutually attracting them- 
selves. The polymer is modeled as a self-avoiding chain 
on a regular (square or simple cubic) lattice with interac- 
tions (euH ! £ffp, epp) = —(1,0,0) between neighboring 
non-bonded monomers. 

This model might be too simple to represent finer de- 
tails of real protein folding , but this is not our concern. 
We use the search for its ground states as a paradig- 
matic example for combinatorial optimization, with a 
large body of existing benchmarks. 

A wide variety of computational strategies have been 
employed to simulate and analyze these models, includ- 
ing conventional (Metropolis) Monte Carlo schemes with 
various types of moves [J, |5|, |(j , chain growth algorithms 
without @ and with re-sampling 00 |r| (see also [ll|). 
genetic algorithms 0, , parallel tempering [l4| and 
generalizations thereof [la, lla] > an 'evolutionary Monte 
Carlo' algorithm [l7j . and others 0]. ^ n addition, Yue 
and Dill [LjJ also devised an exact branch-and-bound al- 
gorithm specific for HP sequences on cubic lattices, which 
gives all low energy states by exact enumeration and typ- 
ically works for N <70 - 80. 

It is the purpose of the present letter to present a 
new variant of the Pruned-Enriched Rosenbluth Method 
(PERM) and to apply it to lattice proteins. PERM is 
a biased chain growth algorithm with re-sampling ( "pop- 
ulation control") and depth-first implementation. It is 
built on the old idea of Rosenbluth and Rosenbluth (RR) 
|2l| to use a biased growth algorithm for polymers, where 
the bias is corrected by means of giving a weight to 
each sample configuration. While the chain grows by 
adding monomers, this weight (which also includes the 
Boltzmann weight if the system is thermal) will fluc- 
tuate. PERM suppresses these fluctuations by "prun- 
ing" configurations with too low weight, and by "enrich- 
ing" the sample with copies of high-weight configurations 
[2fJ. These copies are made while the chain is grow- 
ing, and continue to grow independently of each other. 



PERM can be viewed as a special realization of a "go 
with the winners" strategy |22j and indeed dates back to 
the beginning of the Monte Carlo simulation era, when it 
was called "Russian roulette and splitting" . Among 
statisticians, this approach is also known as sequential 
importance sampling (SIS) with re-sampling 

Pruning and enrichment are done by choosing thresh- 
olds W< and W> depending on the estimate of the par- 
tition sums of n-monomer chains (see below for their ac- 
tual determination). If the current weight W n of an n- 
monomer chain is less than W^, the chain is discarded 
with probability 1/2, otherwise it is kept and its weight 
is doubled. Many alternatives to this simple choice are 
discussed in but we found that more sophisticated 
strategies had little influence on the efficiency, and thus 
we kept the above in the present work. On the contrary, 
we found that different strategies in biasing and, most 
of all, in enrichment had a big effect, and it is here the 
present variant differs from those in 00. There, high- 
weight configurations were simply cloned and the weight 
was uniformly shared between the clones. For relatively 
high temperatures this is very efficient |20j . since each 
clone has so many possibilities to continue that differ- 
ent clones very quickly become independent from each 
other. This is no longer the case for very low temper- 
atures. There we found that clones often evolved in 
the same direction, since one continuation has a much 
higher Boltzmann weight than all others. Thus, cloning 
is no longer efficient in creating configurational diversity, 
which was the main reason why it was introduced. 

The main modification made in the present paper is 
thus that we no longer make identical clones. Rather, 
when we have a configuration with n — 1 monomers, we 
first estimate a predicted weig ht W^ cd for the next step, 
and we count the number kf ree of free sites where the n-th 
monomer can be placed. If kf rcc > 1 and W^ rcd > W> , 
we choose 2 < k < kf ICC different sites among the free 
ones and continue with k configurations which are forced 
to be different. Thus we avoid the loss of diversity which 
limited the success of old PERM. Typically, we used k — 
min{ktr ee ,\Wr d /W>]}. 

When selecting a fc-tuple A = {ai, . . . of mutu- 
ally different continuations ctj with probability pa, the 
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corresponding weights W„ jQl . . . , W n ^ ak are 



IP A 



(1) 



TABLE I: Performances for the 3-d binary (HP-) sequences 
from H. 



where the importance q aj = exp(—(3E n>0lj ) of choice ctj 
is the Boltzmann-Gibbs factor associated with the en- 
ergy B„ )C(j of the newly placed monomer in the poten- 
tial created by all previous monomers. The other terms 
arise from correcting bias and normalization, see |2*ij for 
a more thorough discussion. Choosing 
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would result in usual importance sampling |25j . How- 
ever, instead of q a we use the modified importances 

q a = (*4" e + l/2)g Q in Eq. ©, fc^el being the number 
of free neighbors when the n-th monomer is placed at a. 
This replacement is made since we anticipate that con- 
tinuations with less free neighbours will contribute less 
on the long run than continuations with more free neigh- 
bours. This is similar to "Markovian anticipation" |26j 
within the framework of old PERM, where a bias differ- 
ent from the short-sighted optimal importance sampling 
was found to be preferable. Consequently, the predicted 
weight is W„P rcd = W n - 1 J2 a q a , 

A noteworthy feature of new PERM is that it crosses 
over to complete enumeration when W£ and W> tend 
to zero. In this limit, all possible branches are followed 
and none is pruned as long as its weight is not strictly 
zero. In contrast to this, old PERM would have made 
infinitely many copies of the same configuration. This 
suggests already that we can be more lenient in choosing 
W< and • For the first configuration hitting length 
n we used W„ = and W> = oo, i.e. we neither pruned 
nor branched. For the following configurations we used 
W> = Z n /Z a (c n /c a ) 2 and W< = 0.2 W>. Here, c„ is the 
total number of configurations of length n already created 
during the run, and Z n is the partition sum estimated 
from these configurations. 

In PERM we work at a fixed temperature (no an- 
nealing), and successive "tours" 20] are independent ex- 
cept for the thresholds W^-' y which use partially the 
same partition sum estimates. Results are less sensi- 
tive to the precise choice of temperature than they were 
for old PERM. In general all temperatures in the range 
0.25 < T < 0.35 gave good results for ground state 
search. In the following, when we quote numbers of 
ground state hits or CPU times between such hits, these 
are always independent hits. The actual numbers of (de- 
pendent) hits are much larger. 

We now present our results. Special comparison is 
made with the Core-directed Growth Method (CG) of 
Beutler and Dill [JjJ, the only method we found to be 
still competitive with ours. We emphasize, however, that 
the CG method works only for the HP model and relies 
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34 


40.5 


3.89 


0.23 
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34 


100.2 


1.99 


0.71 
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33 


284.0 


13.45 


6.57 
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32 


74.7 


5.08 


2.55 
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32 


59.2 


6.60 


1.44 
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32 


144.7 


5.37 


3.35 
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31 


26.6 


2.17 


0.46 
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34 


1420.0 


41.41 


10.53 


10 


33 


18.3 


0.47 


0.08 



"Ground state energies Q|. 

'CPU times (minutes) per independent ground state hit, on 167 
MHz Sun ULTRA I work station; from Ref. 0] 
C CPU times (minutes), same machine 
d CPU times (minutes), same machine 



heavily on heuristics, in contrast to our fully blind gen- 
eral purpose approach. 

(a) We first tested the ten 48-mers from 4]. As with 
old PERM, we could reach lowest energy states for all of 
them, but within much shorter CPU times. For all 10 
chains we used the same temperature, exp(l/T) = 18, 
although we could have optimized CPU times by using 
different temperatures for each chain. 

The CPU times for new PERM in Table [Q are typi- 
cally one order of magnitude smaller than those in jllj . 
except for sequence ^9 whose lowest energy was not hit 
in [llj . Since in a SPARC 1 machine was used which 
is slower by a factor w 10 than the 167 MHz Sun UL- 
TRA I used here, this means that our algorithms have 
comparable speeds. We note that introducing a simple 
configurational bias in new PERM [ijj can already give 
a considerable speedup; in this contribution, however, we 
want to concentrate on blind search. 

(b) Next we studied the two 2-d HP-sequences of 
length = 100 of Ref. [5j. They were originally thought 
to have ground states fitting into a 10 x 10 square with 
energies -44 and -46 l S\, but in [j| configurations fitting 
into this square were found with lower energies. More- 
over, when configurations were allowed to have arbitrary 
shape, even lower energies were found 0, 0, 0] . In the 
present work we studied only configurations of the lat- 
ter type. The lowest energies known by now are -48 [l(l ] 
resp. -50 ^5|. The CPU times needed to find them were 
48 min resp. 50h, on machines with w 500 MHz. In con- 
trast, new PERM needed in average 2.6 min resp. 5.8h 
on a 667 MHz DEC Alpha 21264 between any two hits. 

(c) Several 2-d HP-sequences were introduced in [l2| . 
where the authors tried to fold them using a genetic al- 
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gorithm. Except for the shortest chains they were not 
successful, but putative ground states for all of them 
were found in [3, 0, E| • But for the longest of these 
chains (N — 64), the ground state energy E nlin = —42 
was found in |9j only by means of special tricks which 
amount to non-blind search. With blind search, the low- 
est energy reached by PERM was -39. We should stress 
that PERM as used in [9j was blind for all cases except 
this 64-mer (and when it found E = —49 for the second 
jV = 100 chain of la), in contrast to statements to the 
contrary made in yjfl . 

We now found putative ground states for all chains of 
[l2| with blind search. For the 64-mer the average CPU 
time per hit was ca. 30h on the DEC 21264, which seems 
to be roughly comparable to the CPU times needed in 
OE3 but considerably slower than 0. This sequence 
is particularly difficult for any growth algorithm, and the 
fact that we now found it is particularly noteworthy. 

On the other hand, new PERM was much faster than 
for the sequence with N = 60 of [12 It needed « 10 
seconds on the DEC 21264 to hit E min = -36, and w 0.1 
second to hit E = —35. In contrast, E = —36 was never 
hit in 11], while it took 97 minutes to hit E = — 35. 

(d) A 85-mer 2-d HP sequence was given in , where 
it was claimed to have -Emin = — 52. Using a genetic al- 
gorithm, the authors could find only conformations with 
E > —47. In Ref. , using a newly developed evolution- 
ary Monte Carlo (EMC) method, the authors found the 
putative ground state when assuming large parts of its 
known structure as constraints. This amounts of course 
to non-blind search. Without these constraints, the puta- 
tive ground state was not hit in 0] either, although the 
authors claimed their algorithm to be more efficient than 
all previous ones. We easily found states with E = —52, 
but we also found many conformations with E = —53. 
At exp(l/T) = 90 it took ca. 10 min CPU time between 
successive hits on the Sun ULTRA 1. 

(e) Four 3D HP sequences with N = 58, 103, 124, and 
136 were proposed in [2^, ^ as models for actual pro- 
teins or protein fragments. Low energy states for these 
sequences were searched in ^ij using a newly developed 
and supposedly very efficient algorithm. The energies 
reached in 0| were E = —42, —49, —58 and —65, re- 
spectively. We now found lower energy states after only 
few minutes CPU time, for all four chains. For the longer 
ones, the true ground state energies are indeed much 
lower than those found in [l^, see Table ITT1 

Note the very low temperatures needed to fold the very 
longest chains in an optimal time. If we would be in- 
terested in excited states, higher temperatures would be 
better. For instance, to find E — —66 for the 136-mer 
(which is one unit below the lowest energy reached in 
HU), it took just 2.7 seconds/hit on the DEC 21264 when 
using exp(l/T) = 40. 

(f) The only case where we could not find a known 
ground state is a 3-d HP sequence of length 88 given in 



As shown there, it folds into an irregular /3/a-barrel 
with E min = -72. The difficulties of PERM with this se- 
quence are easily understood by looking at the configura- 
tion shown in [111 ]. The nucleus of the hydrophobic core 
is formed by amino acids #36-53. Before its formation, a 
growth algorithm starting at either end has to form very 
unstable and seemingly unnatural structures which are 
stabilized only by this nucleus, a situation similar to the 
64-mer of Ref. 12]. In order to fold also this chain, we 
would have either to start from the middle of the chain 
(as done in 9] for some sequences) or use some other 
heuristics which help the formation of the hydrophobic 
core. Since we wanted our algorithm to be as general and 
"blind" as possible, we did not incorporate such tricks 

A more detailed discussion of our algorithm, the re- 
sults, and comparison with other methods is given else- 
where [25| . A list containing all sequences for which we 
found new lowest energy configurations is given in Ta- 
ble eh 

In the present paper we presented a new version of 
PERM which is a depth-first implementation of the 
' go- with-the- winners' strategy (or sequential importance 
sampling with re-sampling) . The main improvement over 
old PERM is that we now do not make identical clones of 
high weight (partial) configurations, but we branch such 
that each continuation is forced to be different. We do 
not expect this to have much influence for systems at high 
temperatures, but as we showed, it leads to substantial 
improvement at very low temperatures. 

Comparing our results to previous work, we see that 
we found the known lowest energy states in all cases but 
one. Moreover, whenever we could compare with previ- 
ous CPU times, the comparison was favourable for our 
new algorithms, except for the CG method of Beutler 
and Dill 0] . But we should stress that the latter is very 
specific to HP chains, uses strong heuristics regarding the 
formation of a hydrophobic core, and does not give cor- 
rect Boltzmann weights for excited states. All that is not 
true for our method. 

Although our method could be used for a much wider 
range of applications (see |3jJ for applications of PERM) , 
we presented here only results for heteropolymers with 
two types of monomers and the simplest non-trivial in- 
teractions on the square and simple cubic lattices. But 
we applied it also successfully to the HP model on the 
FCC lattice, to off-lattice heteropolymers, and to lattice 
models with more than 2 types of monomers (to be pub- 
lished) . We hope that our results will also foster applica- 
tions to more realistic protein models. We showed only 
results for lowest energy configurations, but we should 
stress that PERM is not only an optimization algorithm. 
It also gives information on the full thermodynamic be- 
haviour. We skipped this here since finding ground states 
is the most difficult problem in general, and sampling ex- 
cited states is easy compared to it. 
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TABLE II: Newly found lowest energy states for binary sequences with interactions e = (eprpr, £hp, epp) = — (1, 0, 0). 









old 


new 




CPU 


N 


d 


Sequence 
example conformation b 


£ min [Ref.] 


Emii\ 
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eT 


time" 


85 


2 


HAPAHiaPeHiaPsHiaPsHiaPsHPaHaPaHaPaHPH 


-52 [12 


-53 


90 


0.03 


58 


3 


PHPH3PH3P2H2PHPH2PH3PHPHPH2P2H3P2HPHP4HP2HP2H2P2HP2H 
ublfl2urfldrfrbrub2lf3lublbrurdfrubdblbufldblfldr2bdfdlu 


-42fl8] 


-44 


30 


0.19 


103 


3 


P2H2PsH2P2H 2 PHP2HP 7 HP3H2PH2P 6 HP2HPHP2 

HP5 H3 P4 H2 PH2 Pb H2 PiHi PHPs H 5 P2 HP2 
ufrbdflfurdfu2rd2 buruf2ulbluld2burdrubrdhbufldblfulf2rd 
bd2b2uflufd3fururd2fu2ru2ldf2urbl2dbdlbulfru2 


-49[18J 


-54[27] 


60 


3.12 


124 
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P3H3PHP4HPSH2P4H2P2H2P4HP4HP2HP2H2P3H2PHPH3P4H3P6 

H2P2H P2HPH P2H P7H P2H3P4H P3H5P4H2PH PH PH PH 

iirnn^niinl t 1 /Tnn nr Tr 1 in * 1 1 t 1 7 11 ri~f ri i TrnriTrniinn^,niiT' r-> nl nrni nTn 

LL 1 ULL 2, l//L/L J till U ^ Lt / / D (-1/ J i t-LL Ltl 1 IL 1 1 L/L-y / U LLU U v U U 1 / j Lit U 1 U$ LLI 2, 

lf3iirdb2d2luflb2rbrfdrfrubulbuf2U2b2dfrbdf2dldf2U2bdrurbulfl 


-58[18] 


-71 


90 


12.3 


136 
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HP 5 HP4HPH2PH 2 P4HPH3P4HPHPH4PiiHP2HP3HPH2P3H 2 P2HP2 
HPHPHPsHPi H 6 P 3 H 2 P 2 H 3 P3 H 2 PH 5 P g HP4HPHP4 
U2b2'rdl2frbdrdlf2lfr2 brblu3fd2rbubd2r2df3dl2ul2blbrdfrurbldrbul2b 
rdrur f2ur $2 ububdlu2 6d 2 blurul2d2 Idr4ubld2hurubu3 br d 2 / 2 U2 Ida Idbubu 


-65 [18] 


-80 


120 


110 



"hours per independent hit on 667 MHz DEC ALPHA 21264 
,, r=right, l=left, f=forward, b=backward, u=up, d=down 
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